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O ' Abstract 

^vj ' This paper proposes a novel chaotic reaction-diffusion model of cellular tu- 

^^ , mor growth and metastasis. The model is based on the multiscale diffusion 

^ ' cancer-invasion model (MDCM) and formulated by introducing strong nonlin- 

'"^ , ear coupling into the MDCM. The new model exhibits temporal chaotic behav- 

C^^ ' ior (which resembles the classical Lorenz strange attractor) and yet retains all 

^N , the characteristics of the MDCM diffusion model. It mathematically describes 

^ both the processes of carcinogenesis and metastasis, as well as the sensitive 

rr\ . dependence of cancer evolution on initial conditions and parameters. On the 

r ) ' basis of this chaotic tumor-growth model, a generic concept of carcinogenesis 

• , and metastasis is formulated. 

o, 

r^ I Keywords: reaction-diffusion tumor growth model, chaotic attractor, sen- 

' I' sitive dependence on initial tumor characteristics 

^— I ■ 1 Introduction 

>: 

1^ ' Cancer is one of the main causes of morbidity and mortality in the world. There 

^SJ . are several different stages in the growth of a tumor before it becomes so large that 

'^ I it causes the patient to die or reduces permanently their quality of life. Developed 

C^^ ' countries are investing large sums of money into cancer research in order to find 

^T^ ■ cures and improve existing treatments. In comparison to molecular biology, cell 

Q ! biology, and drug delivery research, mathematics has so far contributed relatively 

^ I little to the area [l]. 

A number of mathematical models of avascular (solid) tumor growth were re- 
viewed in [2]. These were generally divided into continuum cell population models 
d I described by diffusion partial differential equations (PDEs) of continuum mechanics 

[2113] combined with chemical kinetics, and discrete cell population models described 
by ordinary differential equations (DDEs). 

On the other hand, in many biological systems it is possible to empirically 
demonstrate the presence of attractors that operate starting from different initial 
conditions [5]. Some of these attractors are points, some are closed curves, while 
the others have non-integer, fractal dimension and are termed "strange attractors" 
[6]. It has been proposed that a prerequisite for proper simulating tumor growth by 
computer is to establish whether typical tumor growth patterns are fractal. The frac- 
tal dimension of tumor outlines was empirically determined using the box- counting 
method f7j. In particular, fractal analysis of a breast carcinoma was performed 
using a morphometric method, which is the box-counting method applied to the 
mammogram as well as to the histologic section of a breast carcinoma [8]. 

If tumor growth is chaotic, this could explain the unreliability of treatment and 
prediction of tumor evolution. More importantly, if chaos is established, this could 



X 



be used to adjust strategies for fighting cancer. Treatment could include some form 
of chaos control and/or anti- control]^ 

In this paper, a plausible chaotic diffusion model of tumor growth and metastasis 
is presented. The approach is a combination of theoretical modelling and empirical 
search. 

2 A multiscale diffusion cancer- invasion model 

Recently, a multiscale diffusion cancer-invasion model (MDCM) was presented in 
O [ini [m [m im [m HU [TBI [TT] , which considers cellular and microenvironmental 
factors simultaneously and interactively. The model was classified as hybrid, since 
a continuum deterministic model (based on a system of reaction-diffusion chemo- 
taxis equations) controls the chemical and extracellular matrix (ECM) kinetics and 
a discrete cellular automata-like model (based on a biased random-walk model) 
controls the cell migration and interaction. The interactions of the tumor cells, 
matrix-metalloproteinases (MMs), matrix-degradative enzymes (MDEs) and oxy- 
gen are described by the four coupled rate PDEs: 
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where n denotes the tumor cell density, / is the MM-concentration, m corresponds 
to the MDE-concentration, and c denotes the oxygen concentration. The four vari- 
ables, n, m, /, c, are all functions of the 3-dimensional spatial variable x and time t. 
All equations represent diffusion except ([2]), which shows only temporal evolution 
of the MM-concentration coupled to the MDE-concentration. Z)„ is the tumor cell 
coefficient, Dm is the MDE coefficient and Dc > is the oxygen diffusion coefficient, 
while X, /i. A, 6, a, 7, (3 are positive constants. The other terms respectively denote: 

xV • (nV/)— haptotaxis; 

/xA^— production of MDE by tumor cell; 

Am— decay of MDE; 

(^m/-degradation of MM by MDE; 

ac— natural decay of oxygen; 

7n— oxygen uptake; and 

/9/— production of oxygen by MM. 

Because of its hybrid nature (cells treated as discrete entities and microenviron- 
mental parameters treated as continuous concentrations), the 4-dimensional (4D) 
model ([I|)-((4l) can be directly linked to experimental measurements of those cellular 
and microenvironmental parameters recognized by cancer biologists as important in 



^The chaotic behavior of a system may be artificially weakened or suppressed if it is undesirable. 
This concept is known as control of chaos. The first and the most important method of chaos control 
is the so-called OGY-method, developed by [19]. However, in recent years, a non-traditional concept 
of anti-control of chaos has emerged. Here, the non-chaotic dynamical system is transformed into 
a chaotic one by small controlled perturbation so that useful properties of a chaotic system can be 
utilized m\M\. 



cancer invasion. Furthermore, the fundamental unit of the model is the cell, and 
the complex collective behavior of the tumor emerges as a consequence of inter- 
actions between factors influencing the life cycle and movement of individual cells 

In order to use realistic parameter values, the system of rate equations ([I])-® 
was non-dimensionalised. The resulting 4D scaled system of rate PDEs [llj is given 
by 

— = dnV^n-pV-{n\/f), (5) 

— = --qmf, (6) 
— — = drnV m + Kn — am, (7) 

— = dc'\7'^c + lyf — Lon — (j)c. (8) 

In [TT] the values of the non-dimensional parameters were given as: 

dn = 0.0005, dm = 0.0005, 4 = 0.5, p = 0.01, rj = 50, 

K = l, o- = 0, z^ = 0.5, u; = 0.57, = 0.025. (9) 

The 4D hybrid PDE-model ([I])-([H) can be seen as a special case of a general 
multi-phase tumor growth PDE ([2] equation (12)), 

^ + V • (V,$,) = V • (A^^) + K{^^, Ci) - pi{^^, Q), (10) 

where for phase z, $i is the volume fraction (X]j '^i = 1)) ^i is the velocity, Di is 
the random motility or diffusion, Ai($j,Cj) is the chemical and phase dependent 
production, and Pi{^i, Ci) is the chemical and phase dependent degradation/death. 
The multi-phase model (jlOp has been derived from the basic conservation equations 
for the different chemical species. 
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+ V • N, = Pi, 



where Ci are the concentrations of the chemical species, subindex a for oxygen, b for 
glucose, c for lactate ion, d for carbon dioxide, e for bicarbonate ion, / for chloride 
ion, and g for hydrogen ion concentration; Nj is the flux of each of the chemical 
species inside the tumor spheroid; and Pi is the net rate of consumption/production 
of the chemical species both by tumor cells and due to the chemical reactions with 
other species. 

In this paper we will search for a temporal 3D cancer chaotic attractor^ 'buried' 
within the 4D hybrid spatio-temporal model (dJ-dH). 



^An attractor is defined as the smallest set which cannot be itself decomposed into two or more 
attractors with distinct basins of attraction (this restriction is necessary since a dynamical system 
may have multiple attractors, each with its own basin of attraction). A chaotic, or strange attractor 
is an attractor that has zero measure in the embedding phase-space and has fractal dimension. 
Trajectories within a strange attractor appear to skip around randomly |20II21) . In particular, the 
celebrated Lorenz attractor is given by the set of ODEs [221 123] 



x = a{y~x), y — bx — y ~ xz, z — xy — cz, (11) 

where x, y and z are dynamical variables, constituting the 3D phase-space of the Lorenz system; 
and a, b and c are the parameters of the system. Originally, Lorenz used this model to describe the 



3 A chaotic multi-scale cancer- invasion model 

From the non-dimensional spatio-temporal AC model ©-([H]), discretization was 
performed by neglecting all the spatial derivatives resulting in the following simple 
4D temporal dynamical system. 

n = 0, (12) 

/ = -Tjmf, (13) 

m = nn — am, (14) 

c = vf — ojn — (pc. (15) 



When simulated, the temporal system (J12p - (jl5p with the set of parameters ([9]) 
exhibits a virtually linear temporal behavior with almost no coupling between the 
four concentrations that have very different quantitative values (all phase plots be- 
tween the four concentrations, not shown here, are virtually one-dimensional). To 
see if a modified version of the system (fT2]) - (fT5]l could lead to a chaotic description 
of tumor growth, four new parameters, a, /?, 7, and 6 were introduced. The resulting 
model is 

n = 0, (16) 

/ = a7]{m - /), (17) 

m = finn + 7(7 — c) — m, (18) 

c = I'fm — ujn — 6(j)c. (19) 

The introduction of the parameters (a, /3, 7, 6) was motivated by the fact that 
tumor cell shape (see Figures 6 and 7) represents a visual manifestation of an un- 
derlying balance of forces and chemical reactions p2]. Specifically, the parameters 
represent the following quantities: 

a = tumor cell volume (proliferation/non-proliferation fraction), 

(3 = glucose level, 

7 = number of tumor cells, 

5 = diffusion from the surface (saturation level). 

A tumor is composed of proliferating (P) and quiescent (or non-proliferating) 
(Q) cells. Tumor cells shift from class P to class Q as the tumor grows in size [33]. 
Model dependence on the ratio of proliferation to non-proliferation is introduced via 
the first parameter, a. The discretization of equation ([5]) leads to cell density being 
modelled as a constant in equation (fT2l) . Accordingly, cell density does not play a 
role in the dynamics. In (J16p - ()19p the cell density is re-introduced into the dynamics 
via the cell number, 7. The importance of introducing 7 also appears in connection 



unpredictable behavior of the weather, where x is the rate of convective overturning (convection 
is the process by which heat is transferred by a moving fluid), y is the horizontal temperature 
overturning, and z is the vertical temperature overturning; the parameters are: a = P— proportional 
to the Prandtl number (ratio of the fluid viscosity of a substance to its thermal conductivity, usually 
set at 10), b = _R— proportional to the Rayleigh number (difference in temperature between the top 
and bottom of the system, usually set at 28), and c = K—a number proportional to the physical 
proportions of the region under consideration (width to height ratio of the box which holds the 
system, usually set at 8/3). The Lorenz system (|3ip has the properties: 

1. Symmetry: (x,y, z) -^ (—a;, —y, z) for all values of the parameters, and 

2. The 2— axis {x = y = Q) is invariant (i.e., all trajectories that start on it also end on it). 



with the cychn-dependent kinase (Cdk) inhibitor p27, the level and activity of which 
increase in response to cell density. Levels and activity of Cdk inhibitor p27 also 
increase with differentiation following loss of adhesion to the ECM [29] . 

The ability to estimate the growth pattern of an individual tumor cell type on 
the basis of morphological measurements should have general applicability in cellular 
investigations, cell-growth kinetics, cell transformation and morphogenesis |34j . 

Cell spreading alone is conducive to proliferation and increases in DNA synthesis, 
indicating that cell morphology is a critical determinant of cell function, at least in 
the presence of optimal growth factors and extracellular matrix (ECM) binding [35] . 
In many cells, the changes in morphology can stimulate cell proliferation through 
integrin-mediated signaling, indicating that cell shape may govern how individual 
cells will respond to chemical signals [36]. 




Figure 1: A 3D Lorenz-like chaotic attractor (see |22[|23j ) from the modified tumor 
growth model (|16p - (ll9p . The attractor effectively couples the MM-concentration /, 
the MDE-concentration m, and the oxygen concentration c in a mask-like fashion. 




Figure 2: The m — f phase plot of the 3D attractor. 



Parameters {a, l3,'j,5), introduced in connection with cancer cells morphology 
and dynamics could also influence the very important factor chromatin associated 



with aggressive tumor phenotype and shorter patient survival time. 

c 




Figure 3: The m — c phase plot of the 3D attractor. 




Figure 4: The f — c phase plot of the 3D attractor. 

For computations, the parameters were set to a = 0.06, /3 = 0.05, 7 = 26.5 
and 5 = 40. Small variation of these chosen values would not affect the qualitative 
behavior of the new temporal model (|16p - ()19|) . Simulations of (J16p - (jl9p . using the 
same initial conditions and the same non-dimensional parameters as before, show 
chaotic behavior in the form of Lorenz-like strange attractor in the 3D (f — m — c) 
subspace of the full 4D {n — f — m — c) phase-space (Figures [IH!]) . 

The spatio-temporal system of rate PDEs corresponding to the system in (15)- 



(18) provides the following multi-scale cancer invasion model. 
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The new tumor-growth model (|20p - ()23p retains all the qualities of the original 
AC model dS])-® plus includes the temporal chaotic 'butterfly'-attractor. This 
chaotic behavior may be a more realistic view on the tumor growth, including 
stochastic-like long-term unpredictability and uncontrollability, as well as sensitive 
dependence of a tumor growth on its initial conditions. 

Based on the new tumor growth model (|2U|) - (j23|) . we have re- formulated the fol- 
lowing generic concept of carcinogenesis and metastasis (see scheme on the following 
page). With the model of strongly coupled PDEs, remodelling of extracellular ma- 
trix (ECM) causes a whole process of the movement of invading cells with increased 
haptotaxis, and at the same time decreased enzyme productions level. This can alter 
chromatin structure, which plays an important role in initiating, propagating and 
terminating cellular response to DNA damage [28] . The effect of haptotaxis in the 
process of cells invading could be modelled with travelling-wave (Fisher) equation 

fia. 
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Celular retraining of cancer stem cells 
and/or activation of positive function of 
cyclin-dependent kinase inhibitor p27 [29] 
and/or decreased expression of SATBl ^30, .31j 



The proposed model (j20p - (123[) describes chaotic behavior relevant to the inva- 
sion of cancer cells (see Figures [5] and [6]). As devices for controlling metastasis/chaos 
we suggest the following processes: Cellular retraining of cancer stem cells and/or 
activation of positive function of cyclin-dependent kinase inhibitor p27 and/or de- 
creased expression of SATBl, which is correlated with aggressive tumor phenotype 
in breast cancer and shorter patient survival time [271 ESI EOl EI] . To connect the 
new parameters (a, /3, 7, 6) to the therapeutic regime we note the fact that significant 
factor of any therapy is tumor re-growth during the rest periods between therapy 
applications, which is again dependent on proliferation fraction dynamics a [33] • 

Also, the age of the tumor may translate to different levels of response to a 
standard therapy. As mentioned earlier, the number of cells, 7, is connected with 











Figure 5: Carcinoma intraductale mamae: magnification at 160 x (left) and 420 



X 



(right). 




Figure 6: Carcinoma ductale in situ mamae with magnification at 80 x (left) and 
carcinoma comedo mamae with magnification at 420 x (right). 

the Cdk inhibitor p27, for which levels and activity increase in response with the cell 
density 7, differentiation following loss of adhesion to the ECM. Cdk inhibitor p27 
regulates cell proliferation, cell motility and apoptosis and is the essential element 
for understanding transduction pathways in the regulation of normal and malignant 
cell proliferation as well as it is new hope for therapeutic intervention |29j . 

Introducing the parameter set (a, (3, 7, 6) into the model of tumor cell morphology 
and function could lead to insight into the relationship between these parameters 
and chromatin. By varying these parameters, it may be possible to predict situations 
that result in dangerous levels of chromatin. High levels of chromatin are associated 
with cancer metastasis and some of the most aggressive cancer types [301I3T] . 

The new temporal model (|16p - (|19p (as well as the corresponding spatio-temporal 
tumor-growth model (f20l) " ([23]) ) is not sensitive to variation of the (a, /?, 7, (^)— values, 
but is sensitive to their corresponding places in the equations. 



4 Conclusion 

A plausible chaotic multi-scale cancer-invasion model has been presented. The new 
model was formulated by introducing nonlinear coupling into the existing hybrid 
multiscale Anderson-Chaplain model. The new model describes chaotic behavior, 
as well as sensitive dependence of a tumor evolution on its initial conditions. On 
the basis of this chaotic tumor-growth model, a generic concept of carcinogenesis 
and metastasis was described. Effective cancer control is reflected in progressive 
reduction in cancer mortality |27j. The proposed model suggests a possible solution 
to carcinogenesis and metastasis, by combining mathematical modelling with latest 
medical discoveries. 
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5 Appendix 

5.1 Vector Operators, PDEs and Diffusion Processes 

5.1.1 The Hamilton operator V 

The core vector differential operator V, called 'nabla' or 'del' operator, having prop- 
erties analogous to those of ordinary vectors, invented by Sir William Rowan Hamil- 
ton (1805-1865), is defined in Cartesian (x,y, 2;)-coordinates as 

^ <9 . d . d ^ 
V = T^i + ttJ + tt'^' 

ox ay oz 

where (i,j,k) are the unit vectors of the coordinate axes {x,y,z), respectively. 

Using Hamilton's V— operator, we can define the following three quantities (fun- 
damental for vector calculus and its various applications): 

The gradient. Let (/){x, y, z) define a smootlo scalar field. Then the gradient of (f), 
denoted V^ or grad (p, is a vector field defined by 

, , ^, 5(/>. d(f). d(f) 
grad (f> = V(j)= --1 + --J + -—k. 
ox oy oz 

The gradient grad !?!> is a vector field which points in the direction of the great- 
est rate of increase of the scalar field (j){x,y,z), and whose magnitude is the 
greatest rate of change. 

The divergence. Let V(j;, y, z) = Vii + V2J + VJjk define a smooth vector field. 
Then the divergence of V, denoted V • V or div\', is a scalar field defined by 

,. ^^ ^ ^^ dVi dV2 ^ dVs 

djv V = V • V = — — + — — + — — . 
ox oy oz 



^'smooth' here means difFerentiable at each point (x, y, z) in a certain region of space 
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The divergence div\'{x, y, z) measures the magnitude of a vector field's source 
or sink at a given point {x,y,z). A vector field that has zero divergence 
everywhere is called solenoidal. 

The curl. If \'{x,y,z) is a smooth vector field as above, then the curl or rota- 
tion of V, denoted V x V or rot V, is a vector field defined by the following 
determinant 

i J k 
rotV = \/ X V 



d 


d 


d 


dx 


dy 


dz 


Vl 


V2 


V3 



The curl shows a vector field's rate of rotation, that is, the direction of the 
axis of rotation and the magnitude of the rotation. It can also be described as 
the circulation density. A vector field which has zero curl everywhere is called 
irrotational. 

5.1.2 The Laplacian operator V^ and the fundamental PDEs 

The Laplacian operator is the dot-product of the Hamilton operator by itself, that 
is V^ = V • V (or, divgrad, the divergence of the gradient), defined by 

. 52 92 52 

V = \ \ . 

dx'^ dy"^ dz"^ 

When applied to a smooth scalar field (p{x,y,z), the Laplacian gives 

VV = divgrad ct> = V- (V<^) = ^ + ^ + ^- (24) 

If the expression (f24|) is equal to zero, we get the (elliptic) Laplace PDEjJ 

^ (p = 0, or TT^ + in + ITT = ^■ 

ax^ ay^ az^ 

The Laplace PDE describes various stationary fields: thermostatic, electrostatic, 
magnetostatic, etc. It is always solved subject to certain boundary conditions. 

If the expression (j24|) is proportional (with the constant c) to the first time 
derivative d4>/dt of the scalar field 0, we get the (parabolic) heat PDE, 

2 , del) d'^ct) d'^cj) d'^cj) dcj) 

The heat PDE describes heat conduction, as well as diffusion processes of various 
biophysical natures. It is always solved subject to certain initial and boundary 
conditions. 

If the expression (I24p is equal to the second time derivative d'^4>/dt^ of the scalar 
field (j), we get the (hyperbolic) wave PDE, 

2 _ d'^4> d^ d^(j) d'^4) _ 52</. 

dt"^ ' 9x2 Qy2 Q~/2 Q-fl 

The wave PDE describes wave processes of various biophysical natures. It is always 
solved subject to certain initial and boundary conditions. 



Note that Laplacian can also be applied to the smooth vector field V(a;,j/,z), producing the 
vector Laplace equation, which is important in Maxwell's electrodynamics. 
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5.1.3 General Diffusion PDE 

Density fluctuations in a material undergoing diffusion are described by the diffusion 
PDE, 

^ = \/-{D{cP,r)V^{r,t)), (26) 

where cp^r, t) denotes the density of the diffusing material at location r = (x, y, z) 
and time t and !?(</>, r) is the collective diffusion coefficient for density / at location r. 
If the diffusion coefficient D depends on the density then the equation is nonlinear, 
otherwise it is linear, li D = 1/c is constant, then the equation reduces to the heat 
equation ([25]) . 

The diffusion equation ()26l) can be derived from the continuity PDE, 



— — h V • f = 0, where f is the flux of the diffusing material, (27) 

which states that a change in density in any part of the system is due to inflow and 
outflow of material into and out of that part of the system. The diffusion equation 
(j26p can be obtained easily from the continuity equation (j27p when combined with 
the phenomenological Pick's first law, which assumes that the flux of the diffusing 
material in any part of the system is proportional to the local density gradient, 

{=-D{4>)V<P{r,t). 

The general (scalar) transport PDE, 

— — h V • f (t, r, (p, V(j)) = g(t, r, (p), where g is called the source, (28) 

describes transport phenomena such as heat transfer, mass transfer, fluid dynam- 
ics, etc. All the transfer processes express a certain conservation principle. In this 
respect, any differential equation addresses a certain quantity as its dependent vari- 
able and thus expresses the balance between the phenomena affecting the evolution 
of this quantity. For example, the temperature of a fluid in a heated pipe is affected 
by convection due to the solid-fluid interface, and due to the fluid-fluid interaction. 
Furthermore, temperature is also diffused inside the fluid. For a steady-state prob- 
lem, with the absence of sources, a differential equation governing the temperature 
will express a balance between convection and diffusion. 

If the dependent variable (scalar or vector field) is denoted by (p, the general 
transport PDE (j28p can be rewritten as 

dp(p 



dt 



+ V-{pu(j)) = V-{DV^) + S^ . (29) 



_ ^\ ^ Convection term Diffusion term Source term 

transient term 

The terms in (j29p have the following meaning: 

the transient term, -&■, accounts for the accumulation of (j) in the concerned 
control volume; 

the convection term, accounts for the transport of (j) due to the existence of the 
velocity field (note the velocity u multiplying </>); 

the diffusion term, V-{D\/(p), accounts for the transport of (p due to its gradients; 

the source term, S^f,, accounts for any sources or sinks that either create or destroy 
(p. Any extra terms that cannot be cast into the convection or diffusion terms 
are considered as source terms. 
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5.2 A paradigm of a strange attractor 

Chaos theory, of which turbulence is the most extreme form, started in 1963, when 
Ed Lorenz from MIT took the Navier-Stokes equations from viscous fluid dynamics 
and reduced them into three first-order coupled nonlinear ODEs, to demonstrate the 
idea of sensitive dependence upon initial conditions and associated chaotic behavior. 

5.2.1 Turbulence and chaos theory 

Viscous fluids evolve according to the nonlinear Navier-Stokes PDEs 

^ + u • Vu + Vp/p = i/V^u + f , (30) 

where u = u{x^,t), (i = 1,2,3) is the fluid 3D velocity, p = p{x^,t) is the pressure 
field, p,v are the fluid density and viscosity coefficient, while f = f(x*,t) is the 
nonlinear external energy source. To simplify the problem, we can impose to f the 
so-called Reynolds condition, (f • u) = e, where e is the average rate of energy 
injection. Fluid dynamicists believe that Navier-Stokes equations (|30p accurately 
describe turbulence. 

5.2.2 'Lorenz mask' attractor 

Turbulence is a spatio-temporal chaotic attractor. In general, an attractor is a set 
of system's states (i.e., points in the system's phase-space), invariant under the 
dynamics, towards which neighboring states in a given basin of attraction asymp- 
totically approach in the course of dynamic evolutionO An attractor is defined as 
the smallest unit which cannot be itself decomposed into two or more attractors with 
distinct basins of attraction. This restriction is necessary since a dynamical system 
may have multiple attractors, each with its own basin of attraction (see, e.g., |20j). 

Conservative systems do not have attractors, since the motion is periodic. For 
dissipative dynamical systems, however, volumes shrink exponentially, so attractors 
have volume in nD phase-space. 

In particular, a stable fixed-point surrounded by a dissipative region is an at- 
tractor known as a map sink^ Regular attractors (corresponding to Lyapunov 
exponents) act as limit cycles, in which trajectories circle around a limiting trajec- 
tory which they asymptotically approach, but never reach. The so-called strange 
attractors^ are bounded regions of phase-space (corresponding to positive Lyapunov 
characteristic exponents) having zero measure in the embedding phase-space and 
a fractal dimension. Trajectories within a strange attractor appear to skip around 
randomly. 

In 1963, Ed Lorenz from MIT was trying to improve weather forecasting. Using 
a primitive computer of those days, he discovered the first chaotic attractor. Lorenz 
used three Cartesian variables, {x,y,z), to define atmospheric convection. Changing 
in time, these variables gave him a trajectory in a (Euclidean) 3D-space. From 
all starts, trajectories settle onto a chaotic, or strange attractor. More precisely, 



^A basin of attraction is a set of points in the system's phase-space, such that initial conditions 
chosen in this set dynamically evolve to a particular attractor. 

®A map sink is a stable fixed-point of a map which, in a dissipative dynamical system, is an 
attractor. 

'^A strange attractor is an attracting set that has zero measure in the embedding phase-space 
and has fractal dimension. Trajectories within a strange attractor appear to skip around randomly. 
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Lorenz reduced the Navier-Stokes equations (|3Up for convective Benard fluid flow 
into three first order coupled nonlinear ODEs and demonstrated with these the idea 
of sensitive dependence upon initial conditions and chaos (see [22} [23]). 



Cold Space 




Warm Earth 



Figure 7: Benard cells, showing a typical vortex of a rolling air, with a warm air 
rising in a ring and a cool air descending in the center (left). A simple model of a 
pair of Benard cells given by the celebrated 'Lorenz-butterfly' (or, 'Lorenz-mask') 
strange attractor (right) [25]. 



The celebrated Lorenz ODEs read 



X = a{y — x), 



y = bx — y 



xz, 



xy 



cz, 



(31) 



where x, y and z are dynamical variables, constituting the 3D phase-space of the 
Lorenz system; and a, h and c are the parameters of the system. Originally, Lorenz 
used this model to describe the unpredictable behavior of the weather, where x is the 
rate of convective overturning (convection is the process by which heat is transferred 
by a moving fluid), y is the horizontal temperature overturning, and z is the vertical 
temperature overturning; the parameters are: a = P— proportional to the Prandtl 
number (ratio of the fluid viscosity of a substance to its thermal conductivity, usually 
set at 10), b = i?— proportional to the Rayleigh number (difference in temperature 
between the top and bottom of the system, usually set at 28), and c = K—a number 
proportional to the physical proportions of the region under consideration (width 
to height ratio of the box which holds the system, usually set at 8/3). The Lorenz 
system (fSTI) has the properties: 

1. Symmetry: {x, y, z) -^ {—x, —y, z) for all values of the parameters, and 

2. The z— axis (x = y = 0) is invariant (i.e., all trajectories that start on it also 
end on it). 

Nowadays it is well-known that the Lorenz model is a paradigm for low-dimensional 
chaos in dynamical systems in meteorology, hydrodynamics, laser physics, supercon- 
ductivity, electronics, oil industry, chemical and biological kinetics, etc. 

The 3D phase-portrait of the Lorenz system ([7]) shows the celebrated ^Lorenz 
mask\ a special type of chaotic, strange, or, fractal attractor (see Figure [7]). It 
depicts the famous ^butterfly effect\ (i.e., sensitive dependence on initial conditions, 
that is, a tiny difference in initial conditions is amplified until two outcomes are 
totally different), so that the long term behavior becomes impossible to predict (e.g., 
long term weather forecasting). The Lorenz mask has the following characteristics: 
(i) trajectory does not intersect itself in three dimensions; (ii) trajectory is not 
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periodic or transient; (iii) general form of the shape does not depend on initial 
conditions; and (iv) exact sequence of loops is very sensitive to the initial conditions. 
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